Progetto TPA parallelismo

Authors

Francesco Datres

Prof. Paolo Rech

Logo Unitn

1. Introduzione e obiettivi del progetto

Nel progetto è stata implementata una formula per descrivere la dissipazione del calore da parte di due sorgenti puntiformi a temperatura costante.

dove T_{i,j}(t) rappresenta la temperatura nel punto (i,j) al tempo t, \alpha è il coefficiente di diffusione termica, e \Delta t è il passo temporale.

1.1 Obiettivi specifici

  • Implementare una versione sequenziale in linguaggio C++ della simulazione di diffusione del calore, usata poi come Benchmark.

  • Parallelizzare il codice usando OpenMP e valutando l’hardware su cui si sta operando.

  • Misurare i tempi di esecuzione delle diverse versioni per matrici di dimensione crescente e calcolare:

    • Speedup = tempo sequenziale / tempo parallelo,

    • Efficienza = speedup / numero di thread.

  • Studiare l’effetto del numero di thread, confrontando prestazioni fino a 6 thread fisici e fino a 12 thread logici, per evidenziare eventuali limiti di scalabilità dovuti all’iperthreading.

  • Analizzare l’impatto della dimensione della matrice

1.2 Matrice di dissipazione

Un esempio delle matrici su cui andremo ad operare è riportato nella figura sottostante, dove possiamo vedere le due sorgenti di calore che influenzano l’ambiente circostante. Per rappresentare la matrice di dissipazione è stato utilizzato un codice Python:

Mostra codice Python per matrice di dissipazione 
import numpy as np
import matplotlib.pyplot as plt

# Legge la matrice dal file generato in C++
with open(r"output.txt", "r") as file:
    data = []
    for line in file:
        values = [float(val) for val in line.strip().split()]
        data.append(values)

matrix = np.array(data)

# Mostra la matrice con colori
plt.imshow(matrix, cmap='inferno')
plt.title("Matrice di dissipazione")
plt.colorbar(label="Valore")
plt.xlabel("Colonna")
plt.ylabel("Riga")
plt.show()

Visualizzazione della matrice di dissipazione

2. Architettura Hardware

Per progettare un codice parallelo efficace è fondamentale considerare l’hardware del sistema su cui verrà eseguito. Nel mio caso, il processore ha le seguenti caratteristiche principali

  • 6 core fisici: unità di elaborazione indipendenti, eseguono calcoli simultanei reali.
  • 12 thread logici: grazie all’Hyper-Threading, ogni core gestisce 2 thread, migliorando l’utilizzo delle risorse.
  • Cache L1: 384 KB – la più vicina al core, accesso rapidissimo.
  • Cache L2: 1,5 MB – buffer intermedio tra L1 e L3.
  • Cache L3: 12 MB – condivisa tra tutti i core, migliora la comunicazione.

In generale, è consigliabile usare 6 thread per operazioni fortemente computazionali, sfruttando al massimo i core fisici mentre per carichi computazionali minori si possono utilizzare fino a 12 thread.

3. Codice sequenziale

Per avere un riferimento di partenza e valutare l’efficacia della parallelizzazione, implementiamo innanzitutto la versione sequenziale della simulazione della diffusione del calore. L’algoritmo segue questi passi:

Mostra codice C++ sequenziale 
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <omp.h>    // solo per omp_get_wtime()

using namespace std;

int main(int argc, char const *argv[])
{
    // Controllo argomenti
    if (argc < 3) {
        cerr << "Uso: " << argv[0]
             << " <dimensione_matrice> <epsilon_convergenza>\n";
        return 1;
    }

    // Parametri da riga di comando
    const int base_dim = atoi(argv[1]);
    const double epsilon = atof(argv[2]);
    const int max_steps = 5000;

    double sum_new = 0.0;
    double sum_old = 0.0;

    // Parametri e dimensioni
    const int dimension = base_dim + 2;
    const double alpha = 0.5;
    const double dt = 0.1;

    // Parametri per misurare il tempo
    double finish = 0;
    double start  = 0;
    double elapsed = 0;

    // Allocazione matrici
    double** matrix = new double*[dimension];
    double** next_matrix = new double*[dimension];
    for (int i = 0; i < dimension; ++i) {
        matrix[i] = new double[dimension]();
        next_matrix[i] = new double[dimension]();
    }

    // Inizializzazione
    for (int i = 0; i < dimension; ++i)
        for (int j = 0; j < dimension; ++j)
            matrix[i][j] = next_matrix[i][j] = 0.0;
    
    matrix[5][5]      = 5.0;
    next_matrix[5][5] = 5.0;
    matrix[20][20]      = 3.0;
    next_matrix[20][20]  = 3.0;

    // File tempi
    ofstream time_file("time_seq.csv");
    time_file << "Run,TimeSeconds\n";

    // Inizio cronometraggio
    start = omp_get_wtime();

    // Iterazione con controllo convergenza
    int k;
    for (k = 0; k < max_steps; ++k) {

        for (int i = 1; i < dimension - 1; ++i) {
            for (int j = 1; j < dimension - 1; ++j) {
                next_matrix[i][j] = matrix[i][j]
                    + alpha * dt * (
                        matrix[i+1][j] +
                        matrix[i-1][j] +
                        matrix[i][j+1] +
                        matrix[i][j-1]
                        - 4.0 * matrix[i][j]
                    );
                sum_new += next_matrix[i][j];
                
            }
        }

        // Reimposta le sorgenti
        next_matrix[5][5] = 5.0;
        next_matrix[20][20] = 3.0;

        // Verifica convergenza
        if (abs(sum_new - sum_old) < epsilon && k > 100)
        {
            break;
        }
        sum_old = sum_new;
        sum_new = 0;


        // Scambia le matrici
        swap(matrix, next_matrix);

        // Fine cronometraggio
        finish = omp_get_wtime();
        elapsed = finish - start;
        time_file << k << "," << (elapsed * 1000.0) << "\n";
        start = finish;

    }

    time_file.close();

    // Output finale
    ofstream out("output_seq.txt");
    for (int i = 1; i < dimension-1; ++i) {
        for (int j = 1; j < dimension-1; ++j) {
            out << matrix[i][j] << " ";
        }
        out << "\n";
    }

    // Deallocazione
    for (int i = 0; i < dimension; ++i) {
        delete[] matrix[i];
        delete[] next_matrix[i];
    }
    delete[] matrix;
    delete[] next_matrix;

    // Stampa passi eseguiti
    cout << "Convergenza raggiunta in " << k << " passi.\n";

    return 0;
}

3.1 Struttura codice sequenziale

  1. Allocazione di due matrici:
    Vengono create due matrici quadrate matrix e next_matrix di dimensione (N × N) in precisione double, inizializzate a zero.

    • La dimensione reale della griglia interna è base_dim + 2, dove i bordi (riga 0, riga N−1, colonna 0, colonna N−1) restano costanti.

    • All’interno, le celle da aggiornare sono quelle con indici i = 1..N−2 e j = 1..N−2.

  2. Condizioni iniziali:

    Sono impostate due sorgenti disse di calore in due punti specifici, che vengono mantenute costanti nel corso della simulazione.

  3. Ciclo principale di simulazione:

    L’algoritmo esegue in modo ciclico l’operazione su ogni cella della matrice, fino a un numero massimo di max_steps oppure fino al raggiungimento della convergenza epsilon. A fine ciclo viene eseguito uno swap delle matrici scambiando i puntatori in modo che la matrice matrix diventa la matrice con i valori al tempo t+1 mentre next_matrix viene poi riutilizzata e sovrascritti al passo successivo.

  4. Misurazione dei tempi (per singolo passo)

    Prima del loop viene inizializzata la variabile start che memorizza il tempo all’inizio del ciclo, in seguito ad ogni iterazione del loop viene misurato il tempo in millisecondi impiegato per ogni k iterazione e salvato in un file .CSV.

4. Codice parallelo con OpenMP

Per sfruttare i 6 core fisici (e, se necessario, fino a 12 thread logici) del sistema, è stata parallelizzata la simulazione usando OpenMP.

Mostra codice C++ parallelo 
#include <iostream>
#include <fstream>
#include <omp.h>
#include <vector>
#include <iomanip>  // per stampare in virgola mobile

using namespace std;

int main(int argc, char** argv) {
    if (argc < 3) {
        cerr << "Uso: " << argv[0] << " <dimensione_matrice> <epsilon>\n";
        return 1;
    }
    if (atoi(argv[1])<20)
    {
        cerr << "Dimensione matrince: " << argv[1] << " < 20, matrice non supportata!";
        return 1;
    }
    
    // Parametri da riga di comando
    const int BASE_N = atoi(argv[1]);
    const double EPSILON  = atof(argv[2]);

    // Aggiungo 2 per il bordo di supporto
    int N = BASE_N + 2;

    const double alpha = 0.5;
    const double dt    = 0.1;
 
    // Alloco le matrici 
    double** matrix = new double*[N];
    double** next_matrix = new double*[N];
    for (int i = 0; i < N; ++i) {
        matrix[i] = new double[N];
        next_matrix[i] = new double[N];
    }

    // Inizializzo le matrici
    for (int i = 0; i < N; ++i) {
        for (int j = 0; j < N; j++){
            matrix[i][j] = 0.0;
            next_matrix[i][j] = 0.0;
        }
    }

    // Imposto le sorgenti fisse
    matrix[5][5]         = 5.0;
    next_matrix[5][5]    = 5.0;
    matrix[20][20]       = 3.0;
    next_matrix[20][20]  = 3.0;

    // Imposto variabili per convergenza
    double sum_old  = 0.0;
    double sum_new  = 0.0;
    const int STEPS = 5000; // Limita l'esecuzione nel caso epsilon è eccessivamente piccolo
    bool flag       = 0;

    // Configuro OpenMP
    omp_set_dynamic(0);
    omp_set_num_threads(6);

    // Inizio cronometraggio
    vector<double> times;
    double last_time = omp_get_wtime();

    //int nth = omp_get_max_threads();
    int nth = 6;
    double* sum_buffers = new double[nth];

    // Regione parallela
    #pragma omp parallel
    {
    int tid = omp_get_thread_num();

    // Numero totale di celle da elaborare (escludendo i bordi)
    int internal_N = N - 2;     // lavoro effettivo: righe e colonne da 1 a N-2
    int total_cells = internal_N * internal_N;

    // Celle di base per ogni thread
    int base_cells = total_cells / nth;
    int extra_cells = total_cells % nth;

    // Celle assegnate a questo thread
    // Le extra_cells sono sicuramente minori di nth quindi posso assegnarle in ordine
    int my_cells = base_cells + (tid < extra_cells ? 1 : 0); 

    // Calcolo l'indice di partenza (offset) globale per questo thread
    int offset = tid * base_cells + min(tid, extra_cells);

    // Converto l'indice lineare di partenza in coordinate 2D (i,j)
    int start_i = offset / internal_N + 1;       // +1 per bordi
    int start_j = offset % internal_N + 1;

    // Inizio ciclo
    for (int step = 0; step < STEPS; ++step) {
        double local_sum_new = 0; // privata

        int i = start_i;
        int j = start_j;

        for (int cell = 0; cell < my_cells; ++cell) {
            next_matrix[i][j] = matrix[i][j]
                + alpha * dt * (
                    matrix[i+1][j] +
                    matrix[i-1][j] +
                    matrix[i][j+1] +
                    matrix[i][j-1] -
                    4.0 * matrix[i][j]
                );
            local_sum_new += next_matrix[i][j];

            // Passa alla prossima cella
            j++;
            if (j == N - 1) {
                j = 1;
                i++;
            }
        }

        sum_buffers[tid] = local_sum_new;

        // sincronizzo tutti i thread
        #pragma omp barrier

        // un solo thread resetta le sorgenti e fa lo swap
        #pragma omp single
        {
            sum_new = 0.0;
            for (int t = 0; t < nth; ++t) {
                sum_new += sum_buffers[t];
            }

            if ((abs(sum_new-sum_old)< EPSILON) && step>100 )
            {
                cout<< "Raggiunta convergenza: " << EPSILON << " in " << step << "steps";
                flag = true;
            }

            sum_old = sum_new;

            // Resetto la variabile per il prossimo step
            sum_new = 0.0;
            next_matrix[5][5]   = 5.0;
            next_matrix[20][20] = 3.0;

            // swap dei puntatori
            swap(matrix, next_matrix);
            double now = omp_get_wtime();
            times.push_back(now - last_time);
            last_time = now;
        }
        
        if (flag)
        {break;}
    
        }
    } // fine parallelo

    // Scrivo i tempi di ogni step su file time.txt
    ofstream time_file("time_par.csv");
    time_file << "Run,TimeSeconds\n";
    for (int i = 0; i < times.size(); i++)
    {
        time_file << i << "," << (times[i] * 1000.0) << "\n";
    }
    time_file.close();

    // Salvo il risultato finale su file
    ofstream fout("output.txt");
    for (int i = 1; i < N-1; ++i) {
        for (int j = 1; j < N-1; ++j) {
            fout << fixed << setprecision(3) << matrix[i][j] << (j+1 < N-1 ? ' ' : '\n');
        }
    }
    fout.close();

    // Dealloco
    for (int i = 0; i < N; ++i) {
        delete[] matrix[i];
        delete[] next_matrix[i];
    }
    delete[] matrix;
    delete[] next_matrix;

    return 0;
}

4.1 Ripartizione del carico in griglia 2D

  • La dimensione reale della matrice con bordi è uguale a N = base_dim + 2 ciò permette di lavorare su tutti i punti della matrice compresi appunti i bordi.

  • La parte interna su cui si opera è di dimensione (N−2)×(N−2). Denominato come internal_N = N − 2.

  • Il numero totale di celle interne da aggiornare a ciascun passo è total_cells = internal_N * internal_N.

Usando nth thread, dividiamo uniformemente in:

base_cells  = total_cells / nth;     // numero minimo di celle/punti per thread
extra_cells = total_cells % nth;     // righe/colonne “avanzanti” da distribuire ai primi extra_cells thread
  • A ogni thread tid (da 0 a nth−1) attribuiamo:
my_cells = base_cells + (tid < extra_cells ? 1 : 0);
offset   = tid * base_cells + min(tid, extra_cells);

  • Offset è l’indice lineare (da 0 a total_cells−1) della cella iniziale del blocco.

  • Per convertire offset in coordinate (start_i, start_j), si fa:

start_i = offset / internal_N + 1;   // +1 per “saltare” il bordo riga 0
start_j = offset % internal_N + 1;   // +1 per “saltare” il bordo colonna 0
  • Ogni thread percorre my_cells celle consecutive: incrementa j++, se j == N−1 (fuori bordo destro) allora j = 1 e i++.

In questo modo:

  • Tutte le celle interne sono coperte una volta sola, senza sovrapposizioni.

  • La contiguità spaziale garantisce che ciascun thread acceda a blocchi contigui di memoria, minimizzando le cache‐miss.

Esempio grafico di matrice suddivisa in una griglia per 6 threads, con extra celle

4.2 Dettagli di implementazione OpenMP

1. Configurazione iniziale

omp_set_dynamic(0);
omp_set_num_threads(6);   // 6 core fisici

Disabilitiamo il dynamic adjustment di OpenMP e fissi­amo il numero di thread a 6 .

2. Buffer per la convergenza

Per evitare l’uso di direttive come #pragma omp atomic sulla variabile sum_new creiamo un array di appoggio:

int nth = omp_get_max_threads();
double* sum_buffers = new double[nth];

Ogni thread scrive il proprio local_sum_new in sum_buffers[tid]; poi solo un thread nella sezione single somma tutti i buffer.

3. Regione parallela

#pragma omp parallel
{
  ...
// Ogni thread scrive nel suo buffer
        sum_buffers[tid] = local_sum_new;

        // Barriera per assicurarsi che tutti abbiano scritto su sum_buffers[tid]
        #pragma omp barrier

        // Sezione single: controllo convergenza, reset somma, swap e cronometraggio
        #pragma omp single
        {
          ...
        }

        // Barriera per assicurarsi che il flag sia visibile a tutti
        #pragma omp barrier
        if (flag) {
            break;
        }
    } // fine loop temporale
}     // fine regione parallela

4.3 Verifica correttezza

  • Confronto risultato sequenziale con parallelo

    Il risultato sequenziale è stato confrontato con il risultato del ciclo parallelo. In particolare è stato utilizzato uno script Python per leggere i file output_seq.txt e output_par.txt, confrontare elemento per elemento e segnalare eventuali discrepanze nei valori con una tolleranza di 10^{-12}.

  • Numero di passi di convergenza

    E’ stato controllato che entrambe le versioni a parità di ordine di convergenza epsilon e dimensione della matrice, convergano con lo stesso numero di passi.

Mostra codice Python verifica correttezza: 
import numpy as np

def leggi_matrice(file_path):
    with open(file_path, 'r') as f:
        righe = f.readlines()
        matrice = [list(map(float, r.strip().split())) for r in righe if r.strip()]
    return np.array(matrice)

def confronta_matrici(mat1, mat2, epsilon=1e-6):
    if mat1.shape != mat2.shape:
        print("Le matrici hanno dimensioni diverse.")
        return

    diff = np.abs(mat1 - mat2)
    max_diff = np.max(diff)
    err_pos  = np.where(diff > epsilon)

    if len(err_pos[0]) == 0:
        print("Le matrici sono uguali entro la tolleranza.")
    else:
        print(f" Differenze rilevate in {len(err_pos[0])} celle.")
        print(f" Differenza massima: {max_diff:.0e}")

if __name__ == "__main__":
    mat_seq = leggi_matrice(r"output_seq.txt")
    mat_par = leggi_matrice(r"output_par.txt")
    confronta_matrici(mat_seq, mat_par)

5. Analisi dei tempi e scalabilità

Per valutare l’efficacia della parallelizzazione sono stati eseguiti una serie di Benchmark variando:

  • Dimensione della matrice (N): 128, 256, 512, 1024

  • Numero di thread (T): 1, 2, 4, 6 (core fisici), 8, 12 (thread logici)

Ogni configurazione è stata provata con un ordine di convergenza epsilon = 0.5 per poi salvare i tempi dei singoli passi nei file CSV:

  • times_seq.csv : tempi (ms) per passo della versione sequenziale

  • times_parallel.csv : tempi (ms) per passo della versione parallela

In seguito leggiamo i tempi CSV, analizzando i tempi notiamo come ci sono dei valori che si discostano in modo importante dalla media, per questo escludiamo gli outliers e poi calcoliamo il tempo medio di esecuzione.

5.1 Analisi tempi codice sequenziale

Importiamo i dati CSV dei tempi dell’algoritmo sequenziale, in particolare analizziamo i dati per un’esecuzione con i seguenti parametri:

  • Dimensione della matrice: 1024

  • Ordine di convergenza: 0.5

Come vediamo abbiamo una media del tempo di esecuzione di 8.27 ms.

5.2 Analisi tempi codice sequenziale

Come per l’algoritmo sequenziale, analizziamo i dati per l’esecuzione del ciclo con in codice parallelo con le seguenti caratteristiche:

  • DImensione della matrice: 1024

  • Ordine di convergenza: 0.5

  • Numero threads: 6

In questo caso il tempo di esecuzione è di 2.79 ms.

5.3 Confronto prestazioni

Confrontiamo ora i tempi di esecuzione per 5000 passi per tre algoritmi:

  • sequenziale

  • parallelo 6 core fisici

  • parallelo 12 core virtuali

5.4 Utilizzo CPU

E’ interessante inoltre notare come varia l’utilizzo della CPU durante l’esecuzione dei tre processi:

  • Con 6 thread, l’occupazione dei core fisici è circa del 70% .
  • Con 12 thread, si ha un utilizzo del 100 %.
  • Nella versione sequenziale, la CPU viene utilizzata al 27%.

Ciò evidenzia ulteriormente come la programmazione in parallelo ci aiuta a sfruttare al meglio l’hardware al massimo delle sue potenzialità.

Utilizzo CPU codice parallelo con 12 threads

Utilizzo CPU codice parallelo 6 threads

Utilizzo CPU codice sequenziale

6. Discussione dei risultati e conclusioni

6.1 Speed up ed efficienza

In seguito al confronto dei tempi totali, sono stati misurati per ciascuna dimensione N = {128, 256, 512, 1024, 2048} i seguenti tempi medi (in secondi):

Speedup ed efficienza per diverse dimensioni.
N Speedup 6 th Efficienza 6 th Speedup 12 th Efficienza 12 th
128 1.14 0.19 0.75 0.06
256 2.44 0.41 1.95 0.16
512 2.66 0.44 2.83 0.24
1024 2.89 0.48 3.27 0.27
2048 2.68 0.45 3.20 0.27

Da questi dati sono stati calcolati gli speed up rispetto alla versione sequenziale e le corrispondenti efficienze.

6.1.1 Speedup

Dove T è il numero di threads mentre N è la dimensione della matrice, in seguito è riportata la tabella che confronta i vari speedup per i diversi algoritmi in funzione della dimensione della matrice e numero di threads:

Speedup per diverse dimensioni della matrice
N Speedup (6 thread) Speedup (12 thread) Parallelo 6 vs 12
128 1.14 0.75 0.66
256 2.44 1.95 0.80
512 2.66 2.83 1.06
1024 2.89 3.27 1.13
2048 2.68 3.20 1.19

6.1.2 Efficienza

Dove T è il numero di threads.

Efficienza per diverse dimensioni della matrice
N Efficienza (6 thread) Efficienza (12 thread)
128 0.19 0.13
256 0.41 0.16
512 0.44 0.24
1024 0.48 0.27
2048 0.45 0.27

6.2 Osservazioni

  • Scalabilità con 6 thread

    • All’aumentare di N, lo speedup cresce quasi linearmente fino a circa 2.89 per N = 2048.
    • L’efficienza si assesta tra 0.41 e 0.48, segno che i 6 core fisici vengono sfruttati al meglio a partire da dimensioni medie-grandi (≥ 256).
    • Per N=128 il guadagno è modesto (1.14×) perché l’overhead di sincronizzazione è comparabile al tempo di calcolo complessivo.

  • Scalabilità con 12 thread

    • Lo speedup migliora ulteriormente (fino a 3.27× su N=2048), ma l’efficienza scende sotto 0.27.

    • Il beneficio aggiuntivo da 6 a 12 thread è più evidente su matrici grandi (N >=1024), quando l’attività di calcolo domina nettamente l’overhead.

  • Dimensione minima consigliata

    • A partire da N=256, la versione a 6 thread supera già il doppio della velocità rispetto allo sequenziale, e per N≥512 è conveniente sfruttare pienamente tutti i core fisici.

    • Al di sotto di N=256, l’overhead di creazione dei thread e le barriere riduce drasticamente l’efficacia della parallelizzazione.

6.3 Conclusioni

Nel complesso, il progetto ha raggiunto l’obiettivo di dimostrare come una corretta suddivisione del carico e un uso mirato delle risorse hardware permettano di scalare efficacemente la simulazione di diffusione del calore. In particolare il confronto tra i tempi di esecuzione, gli speedup e le efficienze ottenute con 6 core fisici e 12 thread logici ha evidenziato chiaramente i vantaggi della parallelizzazione: a partire da matrici di dimensione ≥ 256 × 256, i 6 thread fisici garantiscono un miglioramento fino a 3x rispetto alla versione sequenziale, con un efficienza intorno al 40 - 50%. I 12 thread logici portano ad un ulteriore step in termini di velocità con matrici di dimensioni più grandi (N ≥ 1024), ma con un’efficienza inferiore (20-30%). Questo dimostra che il sovraccarico di sincronizzazione e di inizializzazione del parallelismo sono trascurabili con l’aumento della dimensione della matrice, dove prevale l’efficienza di lavoro dovuta ad una suddivisione corretta dei blocchi di lavoro alle cache L1/L2, aumentando sensibilmente la performance.

L’esperienza ha permesso di osservare in modo pratico come bilanciare la complessità di scrittura del codice parallelo con i guadagni reali in tempo di calcolo, fornendo indicazioni precise su quando e come vale la pena impiegare la parallelizzazione in questo tipo di simulazione.